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We provide a theoretical explanation of the results on the intensity distributions and correlation 
functions obtained from a random beam speckle field in nonlinear bulk waveguides reported in the 
recent publication Y. Bromberg et.al., Nat. Photonics 4, 721 (2010). We study both focusing and 
defocusing case and in the limit of small speckle size (short correlated disordered beam) provide 
analytical asymptotes for the intensity probability distributions at the output facet. Additionally we 
provide a simple relation between the speckle sizes at the input and output of a focusing nonlinear 
waveguide. The results are of practical significance for the nonlinear Hanbury Brown and Twiss 
interferometry both in optical waveguides and Bose-Einstein condensates. 
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I. INTRODUCTION 

Extreme events play an important part in many ar- 
eas of physics. In turbulence they determine the non- 
Gaissian statistics of the tails of the probability distribu- 
tion functions for the properties of random flow [1| and 
in the linear theory of random wave localization the mo- 
menta non-selfaveraging quantities, like e.g. wave trans- 
missivity are determined by rare non-typical realizations 
rather than the typical (localized) ones @, 0] . 

In the context of the nonlinear optics the study of the 
extreme events has recently drawn attention in the con- 
text of the optical rogue waves Q - emerging dynam- 
ical objects of very high amplitude and short lifetime. 
Closely related topics are extreme statistics in Raman 
amplification Q and the supercontinuum generation @. 
The appearance of the rogue waves in optical fibers is 
not necessarily a nonlinear effect and can be observed in 
the linear regime as well Q • Finally the statistics of rare 
events (extreme outages) also determines the probability 
of errors in fiber optical communications with distributed 
amplifier spontaneous emission Q. 

Most of the applications above pertain to the field of 
nonlinear fiber optics. Here we study the emergence 
of the high power optical pulses in the context of the 
nonlinear Hanbury Brown and Twiss (HBT) interfer- 
ometry. The goal of this paper is to explain theoret- 
ically recent experimental results reported in Ref.[9j on 
the non-exponential distribution of the intensity distribu- 
tions of disordered optical field propagating in nonlinear 
bulk waveguides. The linear HBT method was first pro- 
posed in 1950s in astronomy as a means of measuring a 
size of a distance light source (e.g. a star) by measur- 
ing the intensity correlation radius of the received light 
[lOj . The latter is the the typical size of an optical speckle 
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(i.e. bright area) observed when multiple waves emitted 
from a thermal source interfere constructively (Tlj . The 
problem allows classical treatment and when the propa- 
gation is linear one can infer that at a distance Z from a 
source of diameter L the size of the speckle S is given by 
a simple relation S = ZX/L, where A is the wavelength 
of emitted light @ . 

Rcf . [9] was the first publication to our knowledge where 
the HBT interferometry in nonlinear light propagation 
and the resulting speckle distribution were studied in 
bulk AlGaAs waveguides. The system is effectively de- 
scribed by the Nonlinear Schrodinger Equation (NLSE) 
which paves the way to application of the nonlinear HBT 
interferometry not only in the field of optics but weakly 
interacting Bose-Einstein condensates as well [12j |. Two 
principal findings of Ref. Q were that the tail of the 
intensity distribution, P(I), of the speckled field is non- 
exponential (unlike in a linear diffraction [ll|) and the 
intensity correlation radius (speckle size) depends on the 
magnitude of the focusing nonlincarity. The first observa- 
tion signifies the fact that the statistics of the optical filed 
after the nonlinear propagation is no longer Gaussian (or, 
equivalently, the amplitude of the field no longer follows 
the Rayleigh distribution). This effect is also known to 
occur in linear systems when a linear wave propagates 
through a disordered media [l3j . Here however we deal 
with a different type of setup where not only is the non- 
linearity present but also disorder is only in the initial 
conditions (incident beam) and not in the media itself. 
For the focusing nonlinearity in an important 2D configu- 
ration of the system (commonly known as 1+1 geometry) 
this phenomenon was correctly attributed by the authors 
of Ref. Q to generation of bright spatial soliton beamlets 
[Til that now play the role of observed speckles. The tails 
of the intensity distribution are determined by extremely 
rare events when an extremely high power (and narrow 
width) soliton is born from a random beam of a finite 
waist and intensity. Our paper seeks to explain theoret- 
ically the profile and shape of the tails of the intensity 
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probability density functions (PDFs) in the 1+1 geome- 
try using the inverse scattering technique (1ST) for the 
NLSE [15| . We first construct a semi-empirical theory of 
a HBT interferometer in the high-power regime where the 
held is dominated by its soliton component. This theory 
is later corroborated by the results obtained analytically 
from the 1ST in the limit of short correlated source field. 
As for the dependence of the intensity correlation radius 
on the magnitude of nonlinearity observed in [9| this is 
of course to be expected in the nonlinear system. Here 
we also supply a simple theoretical result for the aver- 
age speckle size S in the same short-correlated source 
limit, and derive theoretically the linear scaling of its in- 
verse, with the average intensity Jo of the source 
(or, equivalently, with the nonlinear coefficient). We also 
provide theoretical results in the case of de-focusing non- 
linearity where no bright solitons are observed. In all 
cases we perform full numerical simulations using the pa- 
rameters close to those used in the experimental setup of 
Ref.@ to confirm our analytical predictions. 

Note that the propagation of incoherent fields in non- 
linear optics has been studied previously in various con- 
texts (see e.g. Refs. [Hj]). However most of the systems 
considered usually rely on more complicated non-local or 
non-instantaneous models and the quantities calculated 
are usually not particularly relevant for the HBT prob- 
lem. The local 1+1 NLSE is the simplest model where 
one can get both a clear physical understanding of how 
the nonlinearity affects HBT interferometry and obtain 
some analytical results. Similar problem statement oc- 
curs in the context of fiber optics Ref.[17j, where the 
formation of bright solitons from a disordered input was 
employed to illustrate the concept of soliton discrimina- 
tor as a means of nonlinear optical regeneration. How- 
ever in [13] as opposed to the current paper an opposite 
regime was addressed where the emergence of a soliton 
from a disordered output has a relatively small probabil- 
ity and neither the intensity distribution nor the correla- 
tion function of the output were studied. 



II. THE MODEL 

The dynamics of beam propagation along the direction 
of z axis in the presence of diffraction and nonlinearity 
is given by the nonlinear Schrodinger equation (NLSE) 
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where E is the electrical field, x is the spatial transverse 
coordinate, /3o = (27r/A)no is the propagation constant, 
no is the linear refractive index and n 2 is the nonlinear 
coefficient of the medium. Here we will consider both 
attractive (n 2 > 0) and repulsive (n 2 < 0) cases. We 
will also be using dimensionless soliton units £ = x/L, 
where L is some characteristic width, r = z/Lfj where 
Ld = L 2 /3q is the diffraction length (Rayleigh range) and 



u is the dimensionless field E/yI,I= (|ra 2 | (3 Ld/n ) 1 
Then in new dimensionless units we will have 
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When choosing the model for random disordered in- 
put we opt for the form which mimics as close as pos- 
sible the experimental setup and simulation data from 
Ref.[9(. The physical input is defined by two parame- 
ters: the initial speckle size So and the aperture L (we 
also pick the latter as the normalization width for the 
soliton units above). In our numerical simulations the 
spatial resolution Sx is determined according to Sq to in- 
sure sufficient sampling and the width of computational 
domain L' is set large enough to prevent folding during 
the propagation. The input disordered field is modeled in 
Fourier space by N random low frequency modes, where 
N = L' /Sq. As for the complex amplitudes of these 
random modes, a n , we assume that these form a set of 
of independent identical random variables each having a 
uniformly distributed phase and the amplitude sampled 
from a distribution with the average intensity a 2 . This 
pattern is then filtered by the finite aperture L to mani- 
fest the disordered field at the input facet of the nonlinear 
waveguide. As the number of modes is large the central 
limit theorem applies so that the field above, Eq(x), can 
be considered Gaussian with zero mean and the correla- 
tion function: 



(E (x) E*(x')) = a 



2 sin [tt(x - x')N/L'] 
sin [w(x — x')/L'] 



(3) 



The field is of course non-zero only in the aperture 
window [— L/2, L/2] so the formula above applies only 
to this region. The averaged initial intensity is then 
Io = (\E (x)\ 2 ) = a 2 N. Also when Sq is the smallest 
length scale in the problem we will be using the delta- 
correlated approximation when the r.h.s. of Eq.Q is sub- 
stituted with 2D5(x - x') where D = L'a 2 /2 = S I /2. 
This makes further analytical treatment possible in some 
cases and we will often use it in the paper. 

As for the intensity distribution, P(I), since the initial 
field is Gaussian the statistics of the propagated field in 
the linear case (ri2 = 0) will remain Gaussian. As the 
intensity is the modulus squared of the complex Gaus- 
sian variate its normalized value 1/(1} has an exponen- 
tial distribution, the fact which is well known in the linear 
theory of speckle spectra [llj . Our main task in the sub- 
sequent sections will be to determine the modified shapes 
of the intensity distribution P(I) in the presence of non- 
linearity. Of particular interest is the high intensity tail 
of this distribution determined by rare fluctuations lead- 
ing to the field bursts. 



III. THE PHENOMENOLOGICAL APPROACH 

Let us consider the case of focusing (n 2 > 0) NLSE. 
Then it is known that given enough initial power an arbi- 
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trary initial condition Eq(x) evolves into the combination 
of hyperbolic secant constituents (each corresponding to 
an individual bright soliton) and quasilinear radiation 
15]. Here we will adopt a phcnomenological description 
of the intensity distribution and the correlation proper- 
ties of the output field based on the prescribed form of 
the solution as the sum of statistically independent soli- 
ton pulse shapes with prescribed statistical properties. 

A single soliton solution of @ in soliton units is given 
by: 



soliton, 2 r/, is then given by 



u s (t,£) = 2i] sech [4r] v t + 2rj (f-£o)] 
x exp \-i(2v£+ 2t(v 2 - rj 2 ) - 



(4) 



where parameters 77 and v are related to the soliton's 
amplitude, A, and 'velocity' (i.e. the angle of incidence, 
9) , while the parameters £0 an< 4 4> are the soliton's initial 
position and global phase. The total power of an individ- 
ual soliton (or rather its transversal part in the ir-plane) 
P = J I E 1 2 dx is simply proportional to its amplitude: 
P = 4 V IL = 4 V n /\nML. 

The justification for this approach is presented in the 
following chapter llVl where more rigorous 1ST based anal- 
ysis is performed. It is these soliton constituents that 
contribute to the tails of intensity distribution P(I) as 
the linear radiation quickly disperses away from the aper- 
ture. Let us now assume the regime where the density 
of solitons is not too high so that the average minimal 
distance between the solitons is larger than the average 
width of a soliton - a regime that can be called an asymp- 
totically free regime. Then each soliton contributes inde- 
pendently into the intensity distribution (the interference 
effects are neglected) and we may consider a contribution 
from each individual soliton separately. For a single soli- 
ton solution of the NLSE ((TJ) in the real world units the 
intensity at a given point x is given by: 



I(x;rj,v,x ,T) 
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where 2r/ is the amplitude of the soliton in the dimension- 
less soliton units, 2v is its velocity and Xq is the intial 
position of the soliton center (in /xm). As the input is 
random the amplitude, velocity and the initial position 
of the soliton are also random variables. As for their joint 
probability density function P(v, xq, rf) we shall make an 
assumption which is supported both by numerical sim- 
ulations and the following Zakharov-Shabat eigenvalue 
analysis (see below). Namely we assume that soliton ve- 
locity parameter v is independent from the other vari- 
ables and its distribution is uniform over a symmetric 
interval [—Aw/2, Aw/2]. This immediately means that 
for a soliton emitted from the origin the position shift 
due to the fluctuating velocity is also a uniformly dis- 
tributed random variable in the interval [—A/2, A/2], 
where A = 2/avtL and r = Z/Ld is the propagation 
distance in soliton units. 

The conditional probability density of having the value 
of intensity in the vicinity of / given the amplitude of the 
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where I(x; rj, v, xq, t) is given by Eq.([S|) and the angular 
brackets denote additional averaging with respect to the 
marginal PDF of the initial positions P(xo\i~i). 

In order to perform the averaging analytically we will 
only consider the high intensity tails of the PDF when the 
typical soliton width, L/2rj is much more narrow than the 
width of position distribution, A. Additionally we will 
assume that the propagation distance r is large enough 
so A >> L, and therefore the fluctuations in the soli- 
ton initial position, xq , are negligible when compared to 
these due to the random velocity 2v. After all these as- 
sumptions the result of integration ((BJ) can be presented 
as 



P(I\ m x) = j 
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where rj ;» L/2A and I 3> I677 2 / exp[— 8t]Avt]. 

If the number of produced solitons is n > 1 then it is 
relatively straightforward to derive the average minimal 
distance between the neighbouring solitons which is given 
by A/(n 2 + l). As we have assumed that the width of each 
soliton is much more narrow than the average minimal 
inter-soliton distance this implies that the condition for 
the amplitude is in fact 77 > L(n 2 + 1)/2A. 

Finally to get the marginal intensity distribution P(I) 
we need to average Eq. Q over all realizations of soliton 
amplitude 277. Assuming that the marginal PDF i^(??) is 
known and is the same for all realizations with different 
soliton numbers the result reads: 
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where I > (7/(2Awr) 2 ) ((n) 2 + l) 2 and factor (n) takes 
into account that for each realization all n solitons con- 
tribute equally and independently into the intensity and 
we neglect the effects of interference and soliton collisions 
(i.e. overlap). 



IV. 1ST BASED ANALYSIS 

Many properties of the evolving solution of NLSE , 
including the number of emerging solitons, their ampli- 
tudes, energies and velocities can be established by means 
of so called Zakharov-Shabat spectral problem (ZSSP) 
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Here the complex initial field u(0, £) plays the role of 
the "potential" while the complex eigenvalues £ can have 
both discrete and continuous values. It is the discrete 
spectrum = v n + i rj n which determines the soliton 
part of the solution and in the case of a single soliton 
solution the parameters v and rj are are exactly the ones 
featured in Eq. (jU). 



Poisson fit, <n>-10 79 




TABLE I. The main parameters of the simulations for the 
focusing case 







The size of the aperture, L, /im 


50 


The thickness of the slab, d, /i m 


1.5 


The size of the computational domain, L' , um 


4096 


The total number of points, M ' 


8192 


The number of points resolving the aperture, M 


100 


The total number of random modes, N 


2048 


The maximum propagation distance, z, /im 


8000 


Linear refraction index, no 


3.3 


Propagation constant, /3o, pirn -1 


12.61 


The nonlinear coefficient, ri2, cm 2 /GW 


1.67 x 10 -4 


Initial correlation length, So, M m 


2 


The diffraction length, Ld, mm 


33.9 


Normalization intensity, 7, GW/cm 2 


0.05 


Window used for collecting histograms, A, am 


1024 



FIG. 1. (Color online) The distribution of the number of 
emerging solitons for low intensity (a) and high intensity (b) 
regimes 
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We will start with the numerical Monte Carlo simu- 
lations for the number of emerging solitons as well as 
their amplitudes and phases. In the Monte-Carlo simu- 
lations we took 4000 runs with the parameters of random 
Zakharov-Shabat potential given in Table |U Those were 
chosen to be close to experimental values of Ref. [9] . Wc 
performed two runs with different values of initial average 
power, P — a 2 N d L. The first one corresponds to peak 
power of P =lkW and the second to P = 5kW and we 
will refer to these as "low power" and "high power" runs 
correspondingly (keeping in mind that these labels corre- 
spond to the power of the initial disordered filed). Other 
runs with different values of the input parameters (like 
the initial power and the correlation length) were also 
performed but their results were qualitatively the same 
as in the "high power" or "low power" runs so to keep 
the presentation simple and illustrative we only report 
the results for these two. Because the input is random 
the number of emerging solitons will fluctuate around the 
mean. In Figll]we present a distribution of the number 
of emerging solitons for the values of parameters given in 
Table U together with the corresponding Poisson fits. 

One can clearly see that the number of soliton does 
approximately follow the Poisson distribution for low in- 
tensity (small number of solitons) but this approxima- 
tion breaks down for high power run when the number 
of solitons is higher. Similar results were reported earlier 
in Ref. [ljj in the context of the nonlinear fiber optics. 

Next, shown in FigJ5](as dashed and dotted lines) are 




FIG. 2. (Color online) The marginal PDFs of the real (a) and 
imaginary (b) parts of the ZSSP eigenvalues. 



the numerical marginal PDFs for the real and imagi- 
nary parts of ZS eigenvalues P v (v) and P v {r]) which are 
just scaled distributions for soliton final position and 
amplitude. A separate Monte-Carlo run (not shown) 
has demonstrated that the numerical joint distribution 
P(y, 77) is indeed separable so the marginal distributions 
are sufficient. One can see that the PDF for the real part 
P v (v) is close to uniform with the width Av = 75.8 which 
supports the assumptions made earlier in section IIIII In- 
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deed, assuming that solitons are created in the area local- 
ized by the relative small aperture size L it follows that 
the uniform distribution of the real parts of the eigen- 
values with support Ai> yields the uniform distribution 
of soliton position at the output facet with the support 
A = 2A V (L/Ld)z — 1788/im for the values of parameters 
given in Table [I] Moreover our numerical data also shows 
that the support of the distribution A does not depend 
on the average power P, but rather solely on the input 
correlation length So, where, as expected, shorter input 
correlation distances yield broader distributions for v. 

In the limit of delta-correlated initial filed the results 
above can be confirmed analytically. Indeed one can for- 
mally define a 2D density of states (eigenvalues) of the 
non-Hermitian ZSSP © as 



P(v, = 7 ~ V ^) S (V ~ Vn)) 



(10) 



where summation is performed over all discrete eigenval- 
ues for each realization and the averaging is performed 
over all possible realization. It is clear that the density of 
states p(v,C) is (up to a normalization factor) the prob- 
ability density of having a level in the vicinity of point 
(£, T)). Therefore if one knows the density of states it is 
possible in principle to determine the desired level dis- 
tribution. In Ref. [l8[ this quantity was obtained ana- 
lytically in the limit of the Stratonovich delta-correlated 
potential when the support of the potential L, is large. 
The result reads: 



1 (n/D) cothfa/D) 
P[V,r >> nD sinh 2 ^) 



1 
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One can show [19] that in the strict mathematical limit 
of the Stratonovich white noise the result above is in- 
applicable but it holds for any physical process with a 
symmetric field correlation function of finite but small 
radius, like e.g. (j3j) when Sq is much less than the aper- 
ture length. One immediately notices that the quantity 
p(v, rj) does not depend on the real part of the eigen- 
value, v. Therefore the total number of states with given 
imaginary part, 7/ is infinite, i.e. the probability density 
function P{v,vj) is not normalizable in the w-direction. 
This is again the consequence of idealized nature of the 
white noise and for systems with finite correlation radius 
all the quantities in question are of course finite. This 
is indeed confirmed by the numerical results discussed 
above - as we can see the marginal PDF P v (v) is almost 
flat, it does not depend on the value of the initial power 
and its support diverges in the white noise limit. Thus 
the analytical result, Eq. (fTTj) in the short-correlated limit 
explains both the separability of the eigenvalue distri- 
bution and the flat marginal distribution P(v) observed 
numerically. One can now immediately derive the ana- 
lytic expression for the normalized marginal probability 



_ 2 (rj/D) cothfo/D) 
vW ~ D sinh 2 (rj/D) 
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For the specific parameters of our numerical simulations 
last formula in Ea.(|ll|) yields D = 0.61 for the low power 
run and D — 3.08 for the high power run. The corre- 
sponding analytical curves (fl"2")l are plotted as solid lines 
in Fig.Qb and one can observe a rather good agreement 
with the numerics. 

Plugging the amplitude PDF (fl"2"|) into the expression 
© one obtains 




P es t(I) 



, x coshz coth(x coshz) — 1 , 

f(x) = 2 / dx 

sinh (x coshz) 



(13) 



For the large values of argument the x 3> 1 we have 
/(x) ~ 8xKi(2x) and we obtain the asymptote P(I) <x 
i~ 3 / 4 exp[-(//£) 2 i r ) 1 / 2 ] as the high intensity tail of the 



distribution. It turns out that Eq. (IT3|) provides a remark- 
ably good approximation of the tails of the intensity PDF 
- see Fig |3] where we compare it with the results of the 
direct Monte Carlo simulations of the NLSE propagation 
(again 4000 realizations were used). One can see that 



Low power, numerics 
Low power, theory 
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FIG. 3. (Color online) The PDF of the field intensity obtained 
from the direct numerical simulations of Eq.([2]| and the one 
reconstructed from the 1ST theory via formula (|13l) . The 
average intensity was (I) — 3.28 x 10 _2 GW/cm 2 (low power) 
and (I) = 13.3 x lCT 2 GW/cm 2 (high power). 

our analytical result works rather well for I > 10(7) = 
0.328GW/cm 2 (low power limit) and 1.33GW/cm 2 (high 
power limit). The discrepancies at low intensities are due 
to the fact that i) at low intensity the contribution of 
the radiation (completely ignored in our semi-analytical 
scheme) becomes non-negligible and ii) the dilute soli- 
ton gas approximation is not valid for very low power 
(and hence very wide) solitons that overlap and interfere 
significantly which violates the assumptions used in de- 
riving Ea. lflU)) . For the overlapping solitons the phase 
interference becomes an important effect diminishing the 
soliton contribution into the intensity which explains why 
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the analytical result overestimates the probability of low- 
intensity events. One can also observe that the numerical 
PDF for both high and lower power values has an inter- 
esting structure with the inflection point. This inflection 
point corresponds to the crossover between the regime 
of well separated high-intensity pulses and that of broad 
interfering low-intensity solitons. Finally we plot a —3/2 
slope line as a reference. In Ref. [9j it was suggested that 
this slope is rather universal which would imply some 
universal power-law tails. Our results show that it is not 
so. While it does work well as the best fit in the region 
around the inflection point in the high-power case it is 
well off the mark in the low-power case. Also for the 
high power case one can clearly see a crossover to the 
exponential tail as predicted by Eq. (fT3l) . 

Following Ref.Q we can also introduce a normalized 
intensity autocorrelation function as 



g{Ax) 



J (I(x) I(x + Ax)) dx 
J (I(x)) (I(x + Ax)) dx 



(14) 



where I(x) = ^(a;)! 2 is the fluctuating intensity of the 
beam. For a linear medium it can be shown that g(0) = 2 
and then it falls to 3(00) = 1 over a characteristic length 
scale - the intensity correlation length S (also called the 
"speckle size" ) . In Fig. @] we plot the correlation func- 
tion g(Ax) obtained from the same numerical run as 
the the other statistics. If we compare the limiting val- 
ues of the intensity correlation function with the results 
of Appendix [A] which were obtained using only soliton 
component of the solution one can see that the limit- 
ing value g(oo) is very close to unity while theory gives 
the value (n 2 )/(n) 2 ~ l/(n) 0.92 (for both high and 
low power runs) which is close. As for the opposite 
limit, one can see from Fig.Q that the limiting values 
0(0) = (I 2 ) /(I) 2 2.85 (low power) and w 18.5 (high 
power) are far less than the predictions of Eq. (|A6[) (where 
the moments of 77 were taken from distribution (|12|Q . 
g(0) « 17 (low power) and g(0) ~ 26.5 (high power). 
This discrepancy can be easily explained when one looks 
at Fig. ^ where it is clear that the first two moments 
of the intensity determining g(0) are formed by the re- 
gion / ~ (/) where the nonsoliton part of the radiation 
is important and also soliton pulses are wide enough to 
cover most of the sampling region. However the ana- 
lytical prediction of the Appendix A for the correlation 
length itself, S = L/2a v = L/2D « 40/Ltm (low power) 
and ~ 8/im (high power) holds remarkably well. If we 
recall the definition of the parameter D we get (in the 
short-correlated limit) a simple relation between the ob- 
served correlation length, S, the initial correlation length, 
So, and average initial intensity, Iq: 



S = 



L 2 I 
So Io 



(15) 



So the observed correlation length is inversely propor- 
tional to the initial correlation length and the initial av- 
erage intensity (or power). The latter fact confirms the 
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FIG. 4. (Color online) The intensity correlation function ob- 
tained from the direct numerical simulations of Eq.@. 



experimental results of Ref. |9J. Also the correlation 
length in this regime does not depend on the propaga- 
tion distance, r, (i.e. is saturated) which together with 
the scaling law <?(0) oc r (see Eq. (|A6[0 coincide with the 
results obtained in Q by qualitative considerations. 



V. DEFOCUSING CASE 

Let us now turn attention to the de-focusing case where 
712 < 0. The experimental results of [9] show qualitatively 
different behavior of both the correlation function and 
the intensity PDF. In particular the correlation strength 
g(0) goes down with the intensity of the initial speckle 
field (unlike in the focusing case) and the intensity distri- 
bution also look markedly different. Here we will again 
employ the inverse scattering method to study the re- 
sulting intensity PDF. Similar problem for the defocus- 
ing case was studied previously in [20j using asymptotic 
far field expansion developed by Manakov [ljl [21] . The 
ZSSP for focusing and defocusing case differ only by a 
sign of the potential in the second equation so that in 
the defocusing case one will have 



■ d<P2 t 



(16) 



If we assume zero boundary conditions at infinity for the 
defocusing NLSE no bright (or dark) solitons are formed 
and the far field is formed solely by dispersive waves. 

It is known that asymptotically at large z the field 
intensity I(x,z) in the real world units is given by the 
following formula [U, HI H3 : 



I(x,z) 



IL 



D 



2nz 



In 



L D x 
2L z 



(17) 
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where a(£) is determined via the particular solution of 
([IS]) subject to the following boundary conditions 15]: 



0(0; C) 



0(1; 



o(C) « 
HO 



,-tC 



(18) 



Here the spectral parameter £ is real and we have as- 
sumed for definiteness that the initial support of the 
pulse is [0, L] in the real-world units. Coefficients <z(£) 
and b(Q are called the first and second Jost coefficients 
respectively and satisfy the condition |a| 2 — \b\ 2 = 1. Us- 
ing the invariant imbedding approach already developed 
for the focusing case (see e.g. [22|) we can introduce the 
functions a(C;£) = <pi(£)e*e, &(£;£) = Mt)^*, for 
which we will have the following system of equations: 



da(C,0 

on 



ib((,0e 2ia u(0, o(C;0) = l 



= -*a(C,e)e- 2i «u*(0, 6(C;0) = 



(19) 



The Jost coefficients are recovered as a(£) = a(£, 1) and 
HO — HO !)• Note that as w(£) may be considered Gaus- 
sian with the the correlation function given by ([3]) (in the 
real world units) the phase factor exp(2i££) can be ab- 
sorbed into the definition of the random field u(£) so that 
the statistics of both Jost coefficients becomes indepen- 
dent on the spectral parameter £. From (|17p it immedi- 
ately follows that asymptotically the values of the field 
intensity become uncorrelated, i.e. the correlation func- 
tion g(Ax) tends to unity across the traverse dimension 
of the beam as long as the distance z is large. Let us 
parameterize \a\ and |6| as coshx and sinhx respectively 
and introduce the phase difference between the two Jost 
coefficients: ip — Arg[a] — Arg[6]. Then for these two real 



quantities one obtains a system of equations: 



= -Im[e-^u(0], X(0)=0 
= 2 coth2 X Re[e- iv u(0], 



(20) 



where the value of the initial phase <p(0) is chosen so that 
the derivative <p'(0) is finite (see e.g. [22|]). In the delta- 
correlated limit we obtain from the system (12"U1) (treated 
in the Stratonovich sense) the following Fokker-Planck 
equation for the joint PDF P(x> f, [23 



dP „ d . , ,„ , „. D d 2 P 
- = -^-[coth(2 X )P] + -— - 



2D coth 2 (2x) 



d 2 P 
(21) 

According to (fl7|) . (|18l) the statistics of the intensity is 
determined by the statistics of the quantity In cosh X at 
£ = 1 so we can integrate out the dependence on the an- 
gular variable ip using the periodic boundary conditions 
and make a substitution P(x\ — Y(x\ smn %X f° r the 
resulting marginal distribution of the random variable x- 
The resulting equation reads: 



dY 
d£ 



(D/2) d 



dY 

smh(2x) 

ox 



sinh2x dx 
F( X ;0) = 5(x)/sinh(2 X ) 



(22) 



This equation is known in the theory of stochastic pro- 
cesses @, H3] and it has the solution in quadratures: 



Y( X ;i) 



l 



irD 3 



-D/2 



X' exp(- X ' 2 /2£>) 
Vcosh(2 X ') - cosh(2 X ) 



dx 
(23) 

Finally for the intensity PDF one has the following ex- 
pression: 



P(I) 



2-kz 

Tl, 



exp 



~2nz 


r 


Y 


TTZ 


I 












1. 




To 


I 



(24) 



An interesting feature of the distribution above is that it 
is self-similar in the propagation distance, z, i.e. the PDF 
of the quantity y = zI/ILp is universal and depends only 
on the disorder level D. The tail of this PDF is Gaussian 

Also as mentioned before in this far field limit the field 
values at the different points are uncorrelated so that 
g(Ax) = 1, for Ax > S . 

To test the analytical result above we again choose 
the model parameters close to those considered in exper- 
iments Q . Namely we choose ethanol with an absorbing 



dye as a de-focusing nonlinear medium at the wavelength 
of 552nm and assume the values of parameters given in 
Table IE 

Although in the real experiments the pulse attenua- 
tion was quite high at the length scale considered we 
can still use these parameters for illustrative purposes. 
In our simulations we assumed a fixed averaged value 
of the input power, P — 6W, which correspond to the 
effective nonlinear length Lnl — (I^IAi-W^o) -1 = 
5.33mm. The simulations were performed for three val- 
ues of the propagation distance z comparable to the non- 
linear length. The results are given in FigJS] 

One can clearly see that depending on the propaga- 
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TABLE II. The main parameters of the simulations for the 
defocusing case. 



Parameter 


Value 


The size of the aperture, _L, /xm 


50 


The thickness of the slab, d, cm 


2 


The size of the computational domain, Z/, mm 


8.192 


The total number of points, M r 


16384 


The number of points resolving the aperture, NL 


100 


The total number of random modes, -/V 


4096 


Linear refraction index, no 


1.3 


Propagation constant, /3o, /urn -1 


15.62 


The nonlinear coefficient, ri2, cm 2 /W 


-2.6 x 10~ 8 


Initial correlation length, So, /^m 


2 


The diffraction length, Ld, mm 


39.04 


Normalization intensity, 7, W/cm 2 


82.00 


Window used for collecting histograms, A, /jm 


512 



Far field asymptote 
Linear regime 
z=3mm 
z=5mm 
z=10mm 




VI. CONCLUSION 

To conclude, we have studied both analytically and nu- 
merically the statistics of the intensity distribution of a 
disordered short-correlated pulse propagating in nonlin- 
ear media under conditions close to those experimentally 
observed in Ref.Q. In the limit of the delta-correlated 
pulse, when the initial correlation length (speckle size) is 
much less than the aperture size we provide analytical ex- 
pressions for the intensity distribution for both focusing 
and defocusing media. It turns out that the power-law 
tails reported in 0] are not universal and represent an 
approximate fit to a transitional area followed by an ex- 
ponential asymptote in the focusing case and Gaussian 
asymptote in the de-focusing case. For the latter a uni- 
versal analytical formula for the intensity PDF is given 
in the regime when the propagation length (the length 
of the beam) is larger than the nonlinear length. Also 
in the focusing case a simple expression is given for the 
intensity correlation width, S (formula (I15p ) which re- 
lates it to the initial correlation width (speckle size), So 
and the average intensity of the source Jo confirming the 
results of Rcf . [9] . This formula supplements the relation 
^(O)^ = 2\z/L obtained in Ref. Q using quantitative 
arguments and thus allows one to estimate not only the 
linear size of the object, L, but also its average intensity 
Jo and a correlation radius So (or rather the product of 
the two) when the system is in a soliton regime (high 
power, high number of speckle beamlets). 

We would like to thank Yaron Silberberg for draw- 
ing our attention to the problem and Ehud Altman for 
stimulating discussions. SD would like to thank the De- 
partment of Physics of Complex Systems at Weizmann 
Institute of Science for its warm hospitality. 



FIG. 5. (Color online) The PDF of the self-similar variable y 

for different values of the propagation distance. Appendix A: The derivation of the soliton intensity 

correlation function 



tion distance z there are two regimes with two types of 
statistics. When the propagation distance is less than 
the nonlinear length, Ljv^, the nonlinear term in Eq.([2]) 
can be neglected, i.e. the propagation is linear and as 
mentioned before the intensity PDF is exponential. It 
easy to check that in terms of the self-similar variable y 
the linear PDF is 

When the distance exceeds the nonlinear length the far- 
field asymptotes and (f25| hold and the shape of the 
pdf P(y) becomes universal (with the Gaussian tails). 
The intermediate values of the propagation distance fill 
the gap between the two limiting distributions (as clearly 
seen in Figf5]). 



For a single soliton with the intensity profile ([5]) and 
the uniform position distribution it is possible to calcu- 
late the intensity correlation function g\ provided that 
the typical soliton amplitude, a v is much greater than 
the ratio 1/(4Awt), i.e. all typical soliton pulse realisa- 
tions (as well as the correlation function itself) have the 
width much less than the size of the region where soliton 
are eventually distributed, A. Then one obtains 

(I 1 (x)) = (Ii(x + *x)) = -^( V )=I 1 

and 

16 J 2 / rr 3 \ 

(/i(x)/i(ar + Aa;) = -r ( (q cothg-1)) 

Avt \ sum q / 
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where |x|, |xo|, Ax -C A and q — 2i]Ax/L. So for the 
one-soliton correlation function we obtain 



<7i(Ax) = 4 Avt 



I V 3 



\ sinh 2 



Q 



(q cothg- 



(Al) 



where the averaging in the r.h.s. is performed over the 
amplitude distribution. The relative strength of intensity 
fluctuations for one soliton is given by 



9i (0) 



A Avt (t] 3 ) 



(A2) 



Next, let us consider a train of n solitons with sta- 
tistically independent parameters and random, uniform 
phases. First we consider realizations where exactly n- 
solitons are created. Then we have 



(7(x)) =n(l 1 (x))=nl 1 



and 



(i(x)I(x + Ax)) = n(n - l)(h (x)) (h (x + Ax)) 

+ n(h(x)h(x + Ax) +n(n- l)\{Ei (x)£*(x + Ax))| 2 

where E\ (x) is a single-soliton field Q in the real world 
units. Performing additional averaging over all realiza- 
tions with different number of solitons as well as over the 
spatial coordinate x (denoted by an overbar) we arrive 
at: 



ff(Ax) 



1 /A N M» 1 

Ax + 7^-7^ 
(n) VW W 

f \(E 1 (x)E*(x + Ax)W 
I 2 



(A3) 



Again if we assume that the typical width of a soliton 
is much less than its positional support A one can show 
that the field correlation term is position independent 
and is given by 

i^^Cx + Ax))!^ 



2{r,) 



sech[y] sechhy + q] 1 



-i{q/&TT] ) y 



dy 



If additionally Ax <C min[4r(?7), Lt 1 / 2 ] (i.e. random 
phase shift Acj> — 2vAx/L can be neglected) the above 
simplifies to: 



^KE^EUx + AxW 



■qq 



(77) sinhq 



(A5) 



The strength of fluctuations is given by 



(I 2 ) 

.9(0) = y 



4 Avt (rf) 
W 



7~T (A6) 



3(n) [7])* \W° W 
When the argument of the correlation function is large, 
i.e. the inequality L/2Ax <C cr^ holds the function being 
averaged in Ea. (|Al|) decays much faster than the PDF 
Pr! which can be used to obtain the following asymptote: 



9i (Ax) 



3C(3) AvtP v (0) 
~8 W~ 



s) ,A7) 



where Ax 3> L/2a n is assumed. For the full correla- 
tion function if one assumes additionally that Ax 3> 
L max[l/2cr ?) , y/r] the field correlation contribution (IA4|) 
contains the pre- factor (L/Ax) 4 multiplied by a highly 
oscillating integral, so its contribution is neglected and 
one can write down 



g(Ax) w g(oo) 



9i (Ax) 



(A8) 



with gi (Ax) given by Eq. (|A7p above. 

One can see that the correlation function reaches its 
asymptotic value 3(00) = (n 2 ) / (n) 2 — l/(n) following a 
power law. If the number of emerging solitons follows a 
Poisson distribution (where the variance is equal to the 
mean) the limiting value of the correlation function is 
g{oo) — 1 as in the linear case. For the correlation ra- 
dius S one gets an estimate S s» L/2a n . The latter for- 
mula has a transparent physical meaning: the correlation 
length S is a typical speckle size in any optical system. In 
the regime considered here bright solitons play the role 
of "nonlinear speckles" so the typical speckle size is the 
width of a typical soliton, which is given by L/2a rj . 
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